/*=============================================================================
*
*   Copyright (C) 2020 All rights reserved.
*
*   Filename:		quadrangle_shape_function.hpp
*
*   Author: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Date: 2020-03-22 22:55
*
*   Last Editors: Wang Zhecheng - wangzhecheng@yeah.net
*
*   Last modified:	2020-08-04 01:45
*
*   Description:
*
=============================================================================*/
#ifndef QUADRANGLE_SHAPE_FUNCTION_H
#define QUADRANGLE_SHAPE_FUNCTION_H

#include "../../types.hpp"

namespace geoxfem {
namespace geometry {

template <size_t dim, size_t degree>
struct quadrangle_shape_function {
  static constexpr _real r_5 = .5;
  template <typename Shape>
  static VectorXR getStdNByLocal(const VectorXR&, const Shape&) {
    throw "not implemented yet";
  }
  template <typename Shape>
  static MatrixXR getStddNdlxByLocal(const VectorXR&, const Shape&) {
    throw "not implemented yet";
  }
  template <typename Shape>
  static VectorXR getStdNByGlobal(const VectorXR&, const Shape&) {
    throw "not implemented yet";
  }
  template <typename Shape>
  static VectorXR global2local(const VectorXR&, const Shape&) {
    throw "not implemented yet";
  }
};
/*quadrangle| dimension: 2D, order: 1{{{*/

template <>
template <typename Shape>
VectorXR quadrangle_shape_function<2, 1>::getStdNByLocal(const VectorXR& lcoord,
                                                         const Shape&) {
  Vector4R N;
  _real    c1 = lcoord(0);
  _real    c2 = lcoord(1);
  N << (r1 - c1) * (r1 - c2) * r_25,  //
    (r1 - c1) * (r1 + c2) * r_25,     //
    (r1 + c1) * (r1 + c2) * r_25,     //
    (r1 + c1) * (r1 - c2) * r_25;
  return N;
}
template <>
template <typename Shape>
MatrixXR
quadrangle_shape_function<2, 1>::getStddNdlxByLocal(const VectorXR& lcoord,
                                                    const Shape&) {
  _real c1 = lcoord(0);
  _real c2 = lcoord(1);

  MatrixXR dNdlx(2, 4);
  dNdlx << -r_25 * (r1 - c2),  //
    -r_25 * (r1 + c2),         //
    r_25 * (r1 + c2),          //
    r_25 * (r1 - c2),          //
    -r_25 * (r1 - c1),         //
    r_25 * (r1 - c1),          //
    r_25 * (r1 + c1),          //
    -r_25 * (r1 + c1);
  return dNdlx;
}
template <>
template <typename Shape>
VectorXR
quadrangle_shape_function<2, 1>::getStdNByGlobal(const VectorXR& gcoord,
                                                 const Shape&    shp) {
  Matrix4R m;
  Vector4R v;
  _real    x1 = shp.at(0).x();
  _real    y1 = shp.at(0).y();
  _real    x2 = shp.at(1).x();
  _real    y2 = shp.at(1).y();
  _real    x3 = shp.at(2).x();
  _real    y3 = shp.at(2).y();
  _real    x4 = shp.at(3).x();
  _real    y4 = shp.at(3).y();
  m << 1, 1, 1, 1,   //
    x1, x2, x3, x4,  //
    y1, y2, y3, y4,  //
    x1 * y1, x2 * y2, x3 * y3, x4 * y4;
  v << 1, gcoord[0], gcoord[1], gcoord[0] * gcoord[1];
  return m.inverse() * v;
}
template <>
template <typename Shape>
VectorXR quadrangle_shape_function<2, 1>::global2local(const VectorXR& gcoord,
                                                       const Shape&    shp) {
  Vector2R lcoord;
  MatrixXR lcoords(2, 4);
  lcoords << -1, -1, 1, 1,  //
    -1, 1, 1, -1;
  lcoord << lcoords * getStdNByGlobal(gcoord, shp);
  return lcoord;
}

/*}}}*/
/*quadrangle| dimension: 2D, order: 2{{{*/

template <>
template <typename Shape>
VectorXR quadrangle_shape_function<2, 2>::getStdNByLocal(const VectorXR& lcoord,
                                                         const Shape&) {
  VectorXR N(8);
  _real    c1   = lcoord(0);
  _real    c2   = lcoord(1);
  _real    etam = .25 * (1. - c2);
  _real    etap = .25 * (1. + c2);
  _real    xim  = .25 * (1. - c1);
  _real    xip  = .25 * (1. + c1);

  /* N << 0.25 * (1 - c1) * (1 - c2) * (-c1 - c2 - 1),  // */
  /*   0.5 * (1 - c1) * (1 - c2 * c2),                  // */
  /*   0.25 * (1 - c1) * (1 + c2) * (-c1 + c2 - 1),     // */
  /*   0.5 * (1 - c1 * c1) * (1 + c2),                  // */
  /*   0.25 * (1 + c1) * (1 + c2) * (c1 + c2 - 1),      // */
  /*   0.5 * (1 + c1) * (1 - c2 * c2),                  // */
  /*   0.25 * (1 + c1) * (1 - c2) * (c1 - c2 - 1),      // */
  /*   0.5 * (1 - c1 * c1) * (1 - c2); */

  N << 4. * etam * xim * (-c1 - c2 - 1.),  //
    32. * etam * xim * etap,               //
    4. * etap * xim * (-c1 + c2 - 1.),     //
    32. * xim * xip * etap,                //
    4. * etap * xip * (c1 + c2 - 1.),      //
    32. * etap * xip * etam,               //
    4. * xip * etam * (c1 - c2 - 1.),      //
    32. * xim * xip * etam;
  return N;
}
template <>
template <typename Shape>
MatrixXR
quadrangle_shape_function<2, 2>::getStddNdlxByLocal(const VectorXR& lcoord,
                                                    const Shape&) {
  _real c1   = lcoord(0);
  _real c2   = lcoord(1);
  _real etam = .25 * (1. - c2);
  _real etap = .25 * (1. + c2);
  _real xim  = .25 * (1. - c1);
  _real xip  = .25 * (1. + c1);

  MatrixXR dNdlx(2, 8);
  dNdlx << etam * (2. * c1 + c2),  //
    -8. * etam * etap,             //
    etap * (2. * c1 - c2),         //
    -4. * etap * c1,               //
    etap * (2. * c1 + c2),         //
    8. * etap * etam,              //
    etam * (2. * c1 - c2),         //
    -4. * etam * c1,               //

    xim * (c1 + 2. * c2),  //
    -4. * xim * c2,        //
    xim * (2. * c2 - c1),  //
    8. * xim * xip,        //
    xip * (c1 + 2. * c2),  //
    -4. * xip * c2,        //
    xip * (2. * c2 - c1),  //
    -8. * xim * xip;
  return dNdlx;
}
template <>
template <typename Shape>
VectorXR
quadrangle_shape_function<2, 2>::getStdNByGlobal(const VectorXR& gcoord,
                                                 const Shape&    shp) {
  MatrixXR m(8, 8);
  VectorXR v(8);
  _real    x1 = shp.at(0).x();
  _real    y1 = shp.at(0).y();
  _real    x2 = shp.at(1).x();
  _real    y2 = shp.at(1).y();
  _real    x3 = shp.at(2).x();
  _real    y3 = shp.at(2).y();
  _real    x4 = shp.at(3).x();
  _real    y4 = shp.at(3).y();
  _real    x5 = shp.at(4).x();
  _real    y5 = shp.at(4).y();
  _real    x6 = shp.at(5).x();
  _real    y6 = shp.at(5).y();
  _real    x7 = shp.at(6).x();
  _real    y7 = shp.at(6).y();
  _real    x8 = shp.at(7).x();
  _real    y8 = shp.at(7).y();

  _real x = gcoord(0);
  _real y = gcoord(1);

  m << 1, 1, 1, 1, 1, 1, 1, 1,                                               //
    x1, x2, x3, x4, x5, x6, x7, x8,                                          //
    y1, y2, y3, y4, y5, y6, y7, y8,                                          //
    x1 * x1, x2 * x2, x3 * x3, x4 * x4, x5 * x5, x6 * x6, x7 * x7, x8 * x8,  //
    x1 * y1, x2 * y2, x3 * y3, x4 * y4, x5 * y5, x6 * y6, x7 * y7, x8 * y8,  //
    y1 * y1, y2 * y2, y3 * y3, y4 * y4, y5 * y5, y6 * y6, y7 * y7, y8 * y8,  //
    x1 * x1 * y1, x2 * x2 * y2, x3 * x3 * y3, x4 * x4 * y4, x5 * x5 * y5,
    x6 * x6 * y6, x7 * x7 * y7, x8 * x8 * y8,  //
    x1 * y1 * y1, x2 * y2 * y2, x3 * y3 * y3, x4 * y4 * y4, x5 * y5 * y5,
    x6 * y6 * y6, x7 * y7 * y7, x8 * y8 * y8;

  v << 1, x, y, x * x, x * y, y * y, x * x * y, x * y * y;

  return m.inverse() * v;
}
template <>
template <typename Shape>
VectorXR quadrangle_shape_function<2, 2>::global2local(const VectorXR& gcoord,
                                                       const Shape&    shp) {
  Vector2R lcoord;
  MatrixXR lcoords(2, 8);
  lcoords << -1, -1, -1, 0, 1, 1, 1, 0,  //
    -1, 0, 1, 1, 1, 0, -1, -1;
  lcoord << lcoords * getStdNByGlobal(gcoord, shp);
  return lcoord;
}

/*}}}*/

}  // namespace geometry
}  // namespace geoxfem

#endif /* QUADRANGLE_SHAPE_FUNCTION_H */
